Week 1 work
Arin Vansomphone
2023-07-17
ASA Dynamic Risk Predictions Workshop notes
- JMbayes website (with examples): https://drizopoulos.github.io/JMbayes2/
Part 1: Introduction
Explicit outcomes: multiple longitudinal responses. time-to-events
Implicit outcomes: missing data, random visit times
Research question with AIDS dataset: how strong is the association between CD4 cell count and the risk of death?
Research question with PBC dataset: how strong is the association between bilirubin and the risk of death, can it be used for predictions of survival probabilities
When to use separate analysis (only concerns either survival/longitudinal): * Does treatment affect survival? * Are average longitudinal evolutions different between males and females?
When to use joint analysis: * Complex effect estimation: how strong is the association between the longitudinal evolution of CD4 cell counts and the hazard of death? * Time-varying covariates vs. endogenous * Handling implicit outcomes: focus on longitudinal outcomes with dropouts and random visit times
Part 2: The joint-model
Framework: * Use model to describe the evolution of the covariate *
Fit longitudinal data into cox model 
Think about fitting the green line into the hazard function
Random effects explains all interdependencies in the data
- Longitudinal outcome is independent of TTE outcome
- Repeated measurements in longitudinal are independent of each other
For cox model in joint, requires full specification of joint distribution
- Use spline as a flexible model for baseline hazard
jm() function assumes that the ordering of the patients is the same, and that the time unit is the same
Great explanation of function at 1:10
endogeneous: inside patients, such as bio markers
extrogeneous: outside patients, such as nurse caring for patient, hospital
Part 3: Functional forms
- Instantaneous markers are not always optimal, sometimes other
functional forms are needed to accurately model
- Time dependent
- Cumulative effects (can be weighted)
- functional form is specified by
functional_formsarea()calculates the cumulative effectslope()calculates the time-dependent slope
Part 4: Dynamic predictions
Dynamic prediction: survival probabilities are dynamically updated as additional longitudinal information is recorded
Fitted model estimates conditional survival probabilities
Model accounts all measurements present and before, but not in the future. Obviously becomes more accurate with more measurements
Predictions at 2:05
Functional forms have an effect on predictions, so choose the optimal one according to a model selection problem
Model can predict using new data, does not get re-run
Accuracy for predictions: * Discrimination: classifying high-risk and
low-risk patients * Interested in events occuring in a
medically-relevant interval * sensitivity, AUC and ROC curve * Use
function tvROC() * Calibration: how well the joint model
can accurately predict future events * Use a calibration plot:
calibration_plot() * Brier score: combining discrimination
and calibration * Lower brier score is better * Calculate using
tvBrier() * Take the square root of the brier score to find
the average probability difference (less than 5% is good) * To measure
predictive ability, validation sets are needed
Practical run-through
Dataset: prothro, Liver Cirrhosis dataset
library("JMbayes2")## Loading required package: survival## Loading required package: nlme## Loading required package: GLMMadaptive## Loading required package: splines# linear mixed effects model
lmeFit <- lme(pro ~ time + time:treat, data = prothro, random = ~ time | id)
# cox model
CoxFit <- coxph(Surv(Time, death) ~ treat, data = prothros)
# joint-model
jointFit <- jm(CoxFit, lmeFit, time_var = "time")
# add functional_forms argument to specify functional form
summary(jointFit)##
## Call:
## jm(Surv_object = CoxFit, Mixed_objects = lmeFit, time_var = "time")
##
## Data Descriptives:
## Number of Groups: 488 Number of events: 292 (59.8%)
## Number of Observations:
## pro: 2968
##
## DIC WAIC LPML
## marginal 28075.50 28058.46 -14029.52
## conditional 34725.74 34439.72 -17563.87
##
## Random-effects covariance matrix:
##
## StdDev Corr
## (Intr) 19.3032 (Intr)
## time 4.0524 0.0201
##
## Survival Outcome:
## Mean StDev 2.5% 97.5% P Rhat
## treatprednisone -0.1727 0.1614 -0.4970 0.1315 0.2838 1.0008
## value(pro) -0.0418 0.0040 -0.0496 -0.0343 0.0000 1.0123
##
## Longitudinal Outcome: pro (family = gaussian, link = identity)
## Mean StDev 2.5% 97.5% P Rhat
## (Intercept) 73.3298 0.9817 71.3660 75.2442 0.0000 0.9999
## time 0.4516 0.4738 -0.4656 1.3774 0.3356 1.0259
## time:treatprednisone -0.1292 0.6579 -1.4263 1.1372 0.8496 1.0099
## sigma 17.2312 0.2560 16.7351 17.7553 0.0000 1.0059
##
## MCMC summary:
## chains: 3
## iterations per chain: 3500
## burn-in per chain: 500
## thinning: 1
## time: 8 sec- Rule of thumb: Rhat needs to be below 1.1. If not, increase the number of iterations for MCMC
Extract the data of patient 155
dataP155 <- prothro[prothro$id == 155, ]
dataP155$Time <- dataP155$death <- NULL # delete variable as we want to predict thisNote: if you want another longitudinal biomarker, need to make another linear model for it
Prediction plot
sfit <- predict(jointFit, newdata = dataP155[1, ], process = "event",
times = seq(0, 10, length.out = 51),
return_newdata = TRUE)
sfit## id pro time treat start stop event death Time pred_CIF
## 809 155 38 1.0e-06 prednisone 0 0.2710546 0 0 1e-06 0.0000000
## 809.1 155 38 2.0e-01 prednisone 0 0.2710546 0 0 1e-06 0.1134616
## 809.2 155 38 4.0e-01 prednisone 0 0.2710546 0 0 1e-06 0.1989687
## 809.3 155 38 6.0e-01 prednisone 0 0.2710546 0 0 1e-06 0.2663289
## 809.4 155 38 8.0e-01 prednisone 0 0.2710546 0 0 1e-06 0.3212240
## 809.5 155 38 1.0e+00 prednisone 0 0.2710546 0 0 1e-06 0.3671941
## 809.6 155 38 1.2e+00 prednisone 0 0.2710546 0 0 1e-06 0.4065749
## 809.7 155 38 1.4e+00 prednisone 0 0.2710546 0 0 1e-06 0.4409817
## 809.8 155 38 1.6e+00 prednisone 0 0.2710546 0 0 1e-06 0.4715728
## 809.9 155 38 1.8e+00 prednisone 0 0.2710546 0 0 1e-06 0.4991285
## 809.10 155 38 2.0e+00 prednisone 0 0.2710546 0 0 1e-06 0.5241815
## 809.11 155 38 2.2e+00 prednisone 0 0.2710546 0 0 1e-06 0.5471486
## 809.12 155 38 2.4e+00 prednisone 0 0.2710546 0 0 1e-06 0.5683628
## 809.13 155 38 2.6e+00 prednisone 0 0.2710546 0 0 1e-06 0.5880936
## 809.14 155 38 2.8e+00 prednisone 0 0.2710546 0 0 1e-06 0.6065617
## 809.15 155 38 3.0e+00 prednisone 0 0.2710546 0 0 1e-06 0.6239497
## 809.16 155 38 3.2e+00 prednisone 0 0.2710546 0 0 1e-06 0.6404002
## 809.17 155 38 3.4e+00 prednisone 0 0.2710546 0 0 1e-06 0.6560065
## 809.18 155 38 3.6e+00 prednisone 0 0.2710546 0 0 1e-06 0.6708431
## 809.19 155 38 3.8e+00 prednisone 0 0.2710546 0 0 1e-06 0.6849745
## 809.20 155 38 4.0e+00 prednisone 0 0.2710546 0 0 1e-06 0.6984572
## 809.21 155 38 4.2e+00 prednisone 0 0.2710546 0 0 1e-06 0.7113425
## 809.22 155 38 4.4e+00 prednisone 0 0.2710546 0 0 1e-06 0.7236772
## 809.23 155 38 4.6e+00 prednisone 0 0.2710546 0 0 1e-06 0.7355044
## 809.24 155 38 4.8e+00 prednisone 0 0.2710546 0 0 1e-06 0.7468619
## 809.25 155 38 5.0e+00 prednisone 0 0.2710546 0 0 1e-06 0.7577816
## 809.26 155 38 5.2e+00 prednisone 0 0.2710546 0 0 1e-06 0.7682922
## 809.27 155 38 5.4e+00 prednisone 0 0.2710546 0 0 1e-06 0.7784205
## 809.28 155 38 5.6e+00 prednisone 0 0.2710546 0 0 1e-06 0.7881921
## 809.29 155 38 5.8e+00 prednisone 0 0.2710546 0 0 1e-06 0.7976317
## 809.30 155 38 6.0e+00 prednisone 0 0.2710546 0 0 1e-06 0.8067629
## 809.31 155 38 6.2e+00 prednisone 0 0.2710546 0 0 1e-06 0.8155952
## 809.32 155 38 6.4e+00 prednisone 0 0.2710546 0 0 1e-06 0.8241113
## 809.33 155 38 6.6e+00 prednisone 0 0.2710546 0 0 1e-06 0.8322932
## 809.34 155 38 6.8e+00 prednisone 0 0.2710546 0 0 1e-06 0.8401270
## 809.35 155 38 7.0e+00 prednisone 0 0.2710546 0 0 1e-06 0.8476036
## 809.36 155 38 7.2e+00 prednisone 0 0.2710546 0 0 1e-06 0.8547191
## 809.37 155 38 7.4e+00 prednisone 0 0.2710546 0 0 1e-06 0.8614750
## 809.38 155 38 7.6e+00 prednisone 0 0.2710546 0 0 1e-06 0.8678779
## 809.39 155 38 7.8e+00 prednisone 0 0.2710546 0 0 1e-06 0.8739400
## 809.40 155 38 8.0e+00 prednisone 0 0.2710546 0 0 1e-06 0.8796751
## 809.41 155 38 8.2e+00 prednisone 0 0.2710546 0 0 1e-06 0.8850968
## 809.42 155 38 8.4e+00 prednisone 0 0.2710546 0 0 1e-06 0.8902190
## 809.43 155 38 8.6e+00 prednisone 0 0.2710546 0 0 1e-06 0.8950558
## 809.44 155 38 8.8e+00 prednisone 0 0.2710546 0 0 1e-06 0.8996219
## 809.45 155 38 9.0e+00 prednisone 0 0.2710546 0 0 1e-06 0.9039324
## 809.46 155 38 9.2e+00 prednisone 0 0.2710546 0 0 1e-06 0.9080037
## 809.47 155 38 9.4e+00 prednisone 0 0.2710546 0 0 1e-06 0.9118521
## 809.48 155 38 9.6e+00 prednisone 0 0.2710546 0 0 1e-06 0.9154917
## 809.49 155 38 9.8e+00 prednisone 0 0.2710546 0 0 1e-06 0.9189352
## 809.50 155 38 1.0e+01 prednisone 0 0.2710546 0 0 1e-06 0.9221942
## low_CIF upp_CIF
## 809 0.00000000 0.0000000
## 809.1 0.03717772 0.3141227
## 809.2 0.06794302 0.5199672
## 809.3 0.09523238 0.6602401
## 809.4 0.12073251 0.7536370
## 809.5 0.14155980 0.8112603
## 809.6 0.15629531 0.8578796
## 809.7 0.16938166 0.8918004
## 809.8 0.18643586 0.9169940
## 809.9 0.20281096 0.9375112
## 809.10 0.21801270 0.9530290
## 809.11 0.23221165 0.9648173
## 809.12 0.24556498 0.9738068
## 809.13 0.25829630 0.9802099
## 809.14 0.27019033 0.9843905
## 809.15 0.27960896 0.9876075
## 809.16 0.28866421 0.9901030
## 809.17 0.29756012 0.9920583
## 809.18 0.31031221 0.9935954
## 809.19 0.32447311 0.9948176
## 809.20 0.33417186 0.9958026
## 809.21 0.34374247 0.9966013
## 809.22 0.35310894 0.9972520
## 809.23 0.36248599 0.9977843
## 809.24 0.37620089 0.9982213
## 809.25 0.39025750 0.9985811
## 809.26 0.40469455 0.9988776
## 809.27 0.41472089 0.9991215
## 809.28 0.42350592 0.9993579
## 809.29 0.43238794 0.9995927
## 809.30 0.43703992 0.9996821
## 809.31 0.44134054 0.9998156
## 809.32 0.44544157 0.9999178
## 809.33 0.44937362 0.9999560
## 809.34 0.45684734 0.9999824
## 809.35 0.46577067 0.9999953
## 809.36 0.47491189 0.9999989
## 809.37 0.48492628 0.9999998
## 809.38 0.49512198 1.0000000
## 809.39 0.50533315 1.0000000
## 809.40 0.51556616 1.0000000
## 809.41 0.52582714 1.0000000
## 809.42 0.53642914 1.0000000
## 809.43 0.54636979 1.0000000
## 809.44 0.55662572 1.0000000
## 809.45 0.56691491 1.0000000
## 809.46 0.57725642 1.0000000
## 809.47 0.58769190 1.0000000
## 809.48 0.59826448 1.0000000
## 809.49 0.60901724 1.0000000
## 809.50 0.61999342 1.0000000Plotting prediction plot
plot(sfit)
Predicting the longitudinal outcome
Lpred <- predict(jointFit, newdata = dataP155[1, ],
times = seq(0, 10, length.out = 51),
return_newdata = TRUE)
plot(Lpred)
Combining the plots
plot(Lpred, sfit)
Dynamic predictions
n <- nrow(dataP155)
for (i in seq_len(n)) {
Spred <- predict(jointFit, newdata = dataP155[1:i, ],
process = "event", times = seq(0, 10, length.out = 51),
return_newdata = TRUE)
Lpred <- predict(jointFit, newdata = dataP155[1:i, ],
times = seq(0, 10, length.out = 51),
return_newdata = TRUE)
plot(Lpred,Spred)
}
ROC at 3 years with 1 year window
roc <- tvROC(jointFit, newdata = prothro, Tstart = 3, Dt = 1)
plot(roc)
AUC at 3 years with 1 year window
tvAUC(roc)##
## Time-dependent AUC for the Joint Model jointFit
##
## Estimated AUC: 0.6819
## At time: 4
## Using information up to time: 3 (229 subjects still at risk)Calibration plot at 3 years with 1 year window
calibration_plot(jointFit, newdata = prothro, Tstart = 3, Dt = 1)
Brier score at 3 years with 1 year window
tvBrier(jointFit, newdata = prothro, Tstart = 3, Dt = 1)##
## Prediction Error for the Joint Model 'jointFit'
##
## Estimated Brier score: 0.1015
## At time: 4
## For the 229 subjects at risk at time 3
## Number of subjects with an event in [3, 4): 27
## Number of subjects with a censored time in [3, 4): 14
## Accounting for censoring using model-based weightsConcluding remarks
- Pay attention to:
- Model longitudinal flexibly using splines
- Consider functional form for joint model
Tuning longitudinal data through splines
Spline overview: https://rpubs.com/alecri/review_longitudinal
Linear spline: divide axis into segments and consider piecewise-linear trends.
Locations where lines are tied together are known as knots.
Data: Study of 162 girls on body fat percentage pre-menarcheal and post-menarcheal
# loading packages
library(labelled) # labeling data
library(rstatix) # summary statistics##
## Attaching package: 'rstatix'## The following object is masked from 'package:stats':
##
## filterlibrary(ggpubr) # convenient summary statistics and plots## Loading required package: ggplot2library(GGally) # advanced plot## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2library(car) # useful for anova/wald test## Loading required package: carDatalibrary(Epi) # easy getting CI for model coef/pred
library(lme4) # linear mixed-effects models## Loading required package: Matrix##
## Attaching package: 'lme4'## The following object is masked from 'package:Epi':
##
## factorize## The following object is masked from 'package:GLMMadaptive':
##
## negative.binomial## The following object is masked from 'package:nlme':
##
## lmListlibrary(lmerTest) # test for linear mixed-effects models##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## steplibrary(emmeans) # marginal means##
## Attaching package: 'emmeans'## The following object is masked from 'package:GGally':
##
## pigslibrary(multcomp) # CI for linear combinations of model coef## Loading required package: mvtnorm## Loading required package: TH.data## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:rstatix':
##
## select## The following object is masked from 'package:JMbayes2':
##
## area## The following object is masked from 'package:GLMMadaptive':
##
## negative.binomial##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geyserlibrary(geepack) # generalized estimating equations
library(ggeffects) # marginal effects, adjusted predictions
library(gt) # nice tables
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.2 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ lubridate 1.9.2 ✔ tibble 3.2.1
## ✔ purrr 1.0.1 ✔ tidyr 1.3.0## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::collapse() masks nlme::collapse()
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks rstatix::filter(), stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ tidyr::pack() masks Matrix::pack()
## ✖ dplyr::recode() masks car::recode()
## ✖ dplyr::select() masks MASS::select(), rstatix::select()
## ✖ purrr::some() masks car::some()
## ✖ tidyr::unpack() masks Matrix::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors# load data
load(url("http://alecri.github.io/downloads/data/bodyfat.RData"))
head(bodyfat)## # A tibble: 6 × 6
## id agec agem time pbf occasion
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 9.32 13.2 -3.87 7.94 1
## 2 1 10.3 13.2 -2.86 15.6 2
## 3 1 11.2 13.2 -1.95 13.5 3
## 4 1 12.2 13.2 -1 23.2 4
## 5 1 13.2 13.2 0.0500 10.5 5
## 6 1 14.2 13.2 1.05 20.5 6# trajectories overtime
ggplot(bodyfat, aes(agec, pbf)) +
geom_line(aes(group = id)) +
geom_smooth(se = FALSE, size = 2) +
geom_vline(xintercept = 12.86, linetype = "dashed") +
labs(x = "Age, years", y = "Percent Body Fat")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
# fitting model
bodyfat <- bodyfat %>%
mutate(timepost = pmax(time, 0))
lin_lspl0 <- lmer(pbf ~ time + timepost + (time + timepost | id), data = bodyfat)
summary(lin_lspl0)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: pbf ~ time + timepost + (time + timepost | id)
## Data: bodyfat
##
## REML criterion at convergence: 6062.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7742 -0.5900 -0.0359 0.5946 3.3798
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 45.942 6.778
## time 1.631 1.277 0.29
## timepost 2.750 1.658 -0.54 -0.83
## Residual 9.473 3.078
## Number of obs: 1049, groups: id, 162
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 21.3614 0.5646 161.5587 37.837 < 2e-16 ***
## time 0.4171 0.1572 108.4615 2.654 0.00915 **
## timepost 2.0471 0.2280 132.6814 8.980 2.32e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) time
## time 0.351
## timepost -0.515 -0.872# slopes pre-menarcheal and post-menarcheal
K <- rbind(
"population mean pre-menarcheal slope" = c(0, 1, 0),
"population mean post-menarcheal slope" = c(0, 1, 1)
)
tidy(glht(lin_lspl0, linfct = K), conf.int = TRUE)## # A tibble: 2 × 7
## contrast null.value estimate conf.low conf.high statistic adj.p.value
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 population mean … 0 0.417 0.0674 0.767 2.65 0.0156
## 2 population mean … 0 2.46 2.20 2.73 20.7 0sid <- c(4, 16)
expand.grid(
time = seq(-6, 4.1, 0.05),
id = sid
) %>%
mutate(timepost = pmax(time, 0)) %>%
bind_cols(
indiv_pred = predict(lin_lspl0, newdata = .),
marg_pred = predict(lin_lspl0, newdata = ., re.form = ~ 0)
) %>%
ggplot(aes(time, indiv_pred, group = id, col = factor(id))) +
geom_line(aes(y = marg_pred), col = "black", lwd = 1.5) +
geom_line(aes(linetype = "BLUP")) + #empirical best linear unbiased prediction
geom_point(data = filter(bodyfat, id %in% sid),
aes(y = pbf, shape = factor(id), col = factor(id))) +
guides(shape = "none") +
labs(x = "Time relative to menarche, years", y = "Percent Body Fat",
col = "Id", linetype = "Prediction")
Source: http://users.stat.umn.edu/~helwig/notes/smooth-spline-notes.html
Data: prestige of Canadian occupations from 1971 along with average income
# load data
library(car)
data(Prestige)
head(Prestige)## education income women prestige census type
## gov.administrators 13.11 12351 11.16 68.8 1113 prof
## general.managers 12.26 25879 4.02 69.1 1130 prof
## accountants 12.77 9271 15.70 63.4 1171 prof
## purchasing.officers 11.42 8865 9.11 56.8 1175 prof
## chemists 14.62 8403 11.68 73.5 2111 prof
## physicists 15.64 11030 5.13 77.6 2113 prof# plot data
plot(Prestige$income, Prestige$prestige,
xlab = "Income", ylab = "Prestige")
# fit model
library(npreg)
mod.ss <- with(Prestige, ss(income, prestige))
mod.ss##
## Call:
## ss(x = income, y = prestige)
##
## Smoothing Parameter spar = 0.4923283 lambda = 0.0002148891
## Equivalent Degrees of Freedom (Df) 3.467512
## Penalized Criterion (RSS) 12118.05
## Generalized Cross-Validation (GCV) 127.3134summary(mod.ss)##
## Call:
## ss(x = income, y = prestige)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.991 -7.727 -2.453 7.586 33.006
##
## Approx. Signif. of Parametric Effects:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 63.08 2.228 28.310 0.000e+00 ***
## x 61.01 8.035 7.593 1.807e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approx. Signif. of Nonparametric Effects:
## Df Sum Sq Mean Sq F value Pr(>F)
## s(x) 1.468 2125 1448 11.78 0.0001596 ***
## Residuals 98.532 12118 123
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.09 on 98.53 degrees of freedom
## Multiple R-squared: 0.5949, Adjusted R-squared: 0.5845
## F-statistic: 57.35 on 2.468 and 98.53 DF, p-value: <2e-16plot(mod.ss, xlab = "Income", ylab = "Prestige")
with(Prestige, points(income, prestige))
Review on survival models
Source: https://www.geeksforgeeks.org/survival-analysis-in-r/
Survival analysis is the prediction of events at a specified time.
There are two methods of survival analysis in R: * Kaplan-Meier * Cox proportional hazard model
Kaplan-Meier
Used for censored data, not based on underlying probability distribution
library(survival)
lung## inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
## 1 3 306 2 74 1 1 90 100 1175 NA
## 2 3 455 2 68 1 0 90 90 1225 15
## 3 3 1010 1 56 1 0 90 90 NA 15
## 4 5 210 2 57 1 1 90 60 1150 11
## 5 1 883 2 60 1 0 100 90 NA 0
## 6 12 1022 1 74 1 1 50 80 513 0
## 7 7 310 2 68 2 2 70 60 384 10
## 8 11 361 2 71 2 2 60 80 538 1
## 9 1 218 2 53 1 1 70 80 825 16
## 10 7 166 2 61 1 2 70 70 271 34
## 11 6 170 2 57 1 1 80 80 1025 27
## 12 16 654 2 68 2 2 70 70 NA 23
## 13 11 728 2 68 2 1 90 90 NA 5
## 14 21 71 2 60 1 NA 60 70 1225 32
## 15 12 567 2 57 1 1 80 70 2600 60
## 16 1 144 2 67 1 1 80 90 NA 15
## 17 22 613 2 70 1 1 90 100 1150 -5
## 18 16 707 2 63 1 2 50 70 1025 22
## 19 1 61 2 56 2 2 60 60 238 10
## 20 21 88 2 57 1 1 90 80 1175 NA
## 21 11 301 2 67 1 1 80 80 1025 17
## 22 6 81 2 49 2 0 100 70 1175 -8
## 23 11 624 2 50 1 1 70 80 NA 16
## 24 15 371 2 58 1 0 90 100 975 13
## 25 12 394 2 72 1 0 90 80 NA 0
## 26 12 520 2 70 2 1 90 80 825 6
## 27 4 574 2 60 1 0 100 100 1025 -13
## 28 13 118 2 70 1 3 60 70 1075 20
## 29 13 390 2 53 1 1 80 70 875 -7
## 30 1 12 2 74 1 2 70 50 305 20
## 31 12 473 2 69 2 1 90 90 1025 -1
## 32 1 26 2 73 1 2 60 70 388 20
## 33 7 533 2 48 1 2 60 80 NA -11
## 34 16 107 2 60 2 2 50 60 925 -15
## 35 12 53 2 61 1 2 70 100 1075 10
## 36 1 122 2 62 2 2 50 50 1025 NA
## 37 22 814 2 65 1 2 70 60 513 28
## 38 15 965 1 66 2 1 70 90 875 4
## 39 1 93 2 74 1 2 50 40 1225 24
## 40 1 731 2 64 2 1 80 100 1175 15
## 41 5 460 2 70 1 1 80 60 975 10
## 42 11 153 2 73 2 2 60 70 1075 11
## 43 10 433 2 59 2 0 90 90 363 27
## 44 12 145 2 60 2 2 70 60 NA NA
## 45 7 583 2 68 1 1 60 70 1025 7
## 46 7 95 2 76 2 2 60 60 625 -24
## 47 1 303 2 74 1 0 90 70 463 30
## 48 3 519 2 63 1 1 80 70 1025 10
## 49 13 643 2 74 1 0 90 90 1425 2
## 50 22 765 2 50 2 1 90 100 1175 4
## 51 3 735 2 72 2 1 90 90 NA 9
## 52 12 189 2 63 1 0 80 70 NA 0
## 53 21 53 2 68 1 0 90 100 1025 0
## 54 1 246 2 58 1 0 100 90 1175 7
## 55 6 689 2 59 1 1 90 80 1300 15
## 56 1 65 2 62 1 0 90 80 725 NA
## 57 5 5 2 65 2 0 100 80 338 5
## 58 22 132 2 57 1 2 70 60 NA 18
## 59 3 687 2 58 2 1 80 80 1225 10
## 60 1 345 2 64 2 1 90 80 1075 -3
## 61 22 444 2 75 2 2 70 70 438 8
## 62 12 223 2 48 1 1 90 80 1300 68
## 63 21 175 2 73 1 1 80 100 1025 NA
## 64 11 60 2 65 2 1 90 80 1025 0
## 65 3 163 2 69 1 1 80 60 1125 0
## 66 3 65 2 68 1 2 70 50 825 8
## 67 16 208 2 67 2 2 70 NA 538 2
## 68 5 821 1 64 2 0 90 70 1025 3
## 69 22 428 2 68 1 0 100 80 1039 0
## 70 6 230 2 67 1 1 80 100 488 23
## 71 13 840 1 63 1 0 90 90 1175 -1
## 72 3 305 2 48 2 1 80 90 538 29
## 73 5 11 2 74 1 2 70 100 1175 0
## 74 2 132 2 40 1 1 80 80 NA 3
## 75 21 226 2 53 2 1 90 80 825 3
## 76 12 426 2 71 2 1 90 90 1075 19
## 77 1 705 2 51 2 0 100 80 1300 0
## 78 6 363 2 56 2 1 80 70 1225 -2
## 79 3 11 2 81 1 0 90 NA 731 15
## 80 1 176 2 73 1 0 90 70 169 30
## 81 4 791 2 59 1 0 100 80 768 5
## 82 13 95 2 55 1 1 70 90 1500 15
## 83 11 196 1 42 1 1 80 80 1425 8
## 84 21 167 2 44 2 1 80 90 588 -1
## 85 16 806 1 44 1 1 80 80 1025 1
## 86 6 284 2 71 1 1 80 90 1100 14
## 87 22 641 2 62 2 1 80 80 1150 1
## 88 21 147 2 61 1 0 100 90 1175 4
## 89 13 740 1 44 2 1 90 80 588 39
## 90 1 163 2 72 1 2 70 70 910 2
## 91 11 655 2 63 1 0 100 90 975 -1
## 92 22 239 2 70 1 1 80 100 NA 23
## 93 5 88 2 66 1 1 90 80 875 8
## 94 10 245 2 57 2 1 80 60 280 14
## 95 1 588 1 69 2 0 100 90 NA 13
## 96 12 30 2 72 1 2 80 60 288 7
## 97 3 179 2 69 1 1 80 80 NA 25
## 98 12 310 2 71 1 1 90 100 NA 0
## 99 11 477 2 64 1 1 90 100 910 0
## 100 3 166 2 70 2 0 90 70 NA 10
## 101 1 559 1 58 2 0 100 100 710 15
## 102 6 450 2 69 2 1 80 90 1175 3
## 103 13 364 2 56 1 1 70 80 NA 4
## 104 6 107 2 63 1 1 90 70 NA 0
## 105 13 177 2 59 1 2 50 NA NA 32
## 106 12 156 2 66 1 1 80 90 875 14
## 107 26 529 1 54 2 1 80 100 975 -3
## 108 1 11 2 67 1 1 90 90 925 NA
## 109 21 429 2 55 1 1 100 80 975 5
## 110 3 351 2 75 2 2 60 50 925 11
## 111 13 15 2 69 1 0 90 70 575 10
## 112 1 181 2 44 1 1 80 90 1175 5
## 113 10 283 2 80 1 1 80 100 1030 6
## 114 3 201 2 75 2 0 90 100 NA 1
## 115 6 524 2 54 2 1 80 100 NA 15
## 116 1 13 2 76 1 2 70 70 413 20
## 117 3 212 2 49 1 2 70 60 675 20
## 118 1 524 2 68 1 2 60 70 1300 30
## 119 16 288 2 66 1 2 70 60 613 24
## 120 15 363 2 80 1 1 80 90 346 11
## 121 22 442 2 75 1 0 90 90 NA 0
## 122 26 199 2 60 2 2 70 80 675 10
## 123 3 550 2 69 2 1 70 80 910 0
## 124 11 54 2 72 1 2 60 60 768 -3
## 125 1 558 2 70 1 0 90 90 1025 17
## 126 22 207 2 66 1 1 80 80 925 20
## 127 7 92 2 50 1 1 80 60 1075 13
## 128 12 60 2 64 1 1 80 90 993 0
## 129 16 551 1 77 2 2 80 60 750 28
## 130 12 543 1 48 2 0 90 60 NA 4
## 131 4 293 2 59 2 1 80 80 925 52
## 132 16 202 2 53 1 1 80 80 NA 20
## 133 6 353 2 47 1 0 100 90 1225 5
## 134 13 511 1 55 2 1 80 70 NA 49
## 135 1 267 2 67 1 0 90 70 313 6
## 136 22 511 1 74 2 2 60 40 96 37
## 137 12 371 2 58 2 1 80 70 NA 0
## 138 13 387 2 56 1 2 80 60 1075 NA
## 139 1 457 2 54 1 1 90 90 975 -5
## 140 5 337 2 56 1 0 100 100 1500 15
## 141 21 201 2 73 2 2 70 60 1225 -16
## 142 3 404 1 74 1 1 80 70 413 38
## 143 26 222 2 76 1 2 70 70 1500 8
## 144 1 62 2 65 2 1 80 90 1075 0
## 145 11 458 1 57 1 1 80 100 513 30
## 146 26 356 1 53 2 1 90 90 NA 2
## 147 16 353 2 71 1 0 100 80 775 2
## 148 16 163 2 54 1 1 90 80 1225 13
## 149 12 31 2 82 1 0 100 90 413 27
## 150 13 340 2 59 2 0 100 90 NA 0
## 151 13 229 2 70 1 1 70 60 1175 -2
## 152 22 444 1 60 1 0 90 100 NA 7
## 153 5 315 1 62 2 0 90 90 NA 0
## 154 16 182 2 53 2 1 80 60 NA 4
## 155 32 156 2 55 1 2 70 30 1025 10
## 156 NA 329 2 69 1 2 70 80 713 20
## 157 26 364 1 68 2 1 90 90 NA 7
## 158 4 291 2 62 1 2 70 60 475 27
## 159 12 179 2 63 1 1 80 70 538 -2
## 160 1 376 1 56 2 1 80 90 825 17
## 161 32 384 1 62 2 0 90 90 588 8
## 162 10 268 2 44 2 1 90 100 2450 2
## 163 11 292 1 69 1 2 60 70 2450 36
## 164 6 142 2 63 1 1 90 80 875 2
## 165 7 413 1 64 1 1 80 70 413 16
## 166 16 266 1 57 2 0 90 90 1075 3
## 167 11 194 2 60 2 1 80 60 NA 33
## 168 21 320 2 46 1 0 100 100 860 4
## 169 6 181 2 61 1 1 90 90 730 0
## 170 12 285 2 65 1 0 100 90 1025 0
## 171 13 301 1 61 1 1 90 100 825 2
## 172 2 348 2 58 2 0 90 80 1225 10
## 173 2 197 2 56 1 1 90 60 768 37
## 174 16 382 1 43 2 0 100 90 338 6
## 175 1 303 1 53 1 1 90 80 1225 12
## 176 13 296 1 59 2 1 80 100 1025 0
## 177 1 180 2 56 1 2 60 80 1225 -2
## 178 13 186 2 55 2 1 80 70 NA NA
## 179 1 145 2 53 2 1 80 90 588 13
## 180 7 269 1 74 2 0 100 100 588 0
## 181 13 300 1 60 1 0 100 100 975 5
## 182 1 284 1 39 1 0 100 90 1225 -5
## 183 16 350 2 66 2 0 90 100 1025 NA
## 184 32 272 1 65 2 1 80 90 NA -1
## 185 12 292 1 51 2 0 90 80 1225 0
## 186 12 332 1 45 2 0 90 100 975 5
## 187 2 285 2 72 2 2 70 90 463 20
## 188 3 259 1 58 1 0 90 80 1300 8
## 189 15 110 2 64 1 1 80 60 1025 12
## 190 22 286 2 53 1 0 90 90 1225 8
## 191 16 270 2 72 1 1 80 90 488 14
## 192 16 81 2 52 1 2 60 70 1075 NA
## 193 12 131 2 50 1 1 90 80 513 NA
## 194 1 225 1 64 1 1 90 80 825 33
## 195 22 269 2 71 1 1 90 90 1300 -2
## 196 12 225 1 70 1 0 100 100 1175 6
## 197 32 243 1 63 2 1 80 90 825 0
## 198 21 279 1 64 1 1 90 90 NA 4
## 199 1 276 1 52 2 0 100 80 975 0
## 200 32 135 2 60 1 1 90 70 1275 0
## 201 15 79 2 64 2 1 90 90 488 37
## 202 22 59 2 73 1 1 60 60 2200 5
## 203 32 240 1 63 2 0 90 100 1025 0
## 204 3 202 1 50 2 0 100 100 635 1
## 205 26 235 1 63 2 0 100 90 413 0
## 206 33 105 2 62 1 2 NA 70 NA NA
## 207 5 224 1 55 2 0 80 90 NA 23
## 208 13 239 2 50 2 2 60 60 1025 -3
## 209 21 237 1 69 1 1 80 70 NA NA
## 210 33 173 1 59 2 1 90 80 NA 10
## 211 1 252 1 60 2 0 100 90 488 -2
## 212 6 221 1 67 1 1 80 70 413 23
## 213 15 185 1 69 1 1 90 70 1075 0
## 214 11 92 1 64 2 2 70 100 NA 31
## 215 11 13 2 65 1 1 80 90 NA 10
## 216 11 222 1 65 1 1 90 70 1025 18
## 217 13 192 1 41 2 1 90 80 NA -10
## 218 21 183 2 76 1 2 80 60 825 7
## 219 11 211 1 70 2 2 70 30 131 3
## 220 2 175 1 57 2 0 80 80 725 11
## 221 22 197 1 67 1 1 80 90 1500 2
## 222 11 203 1 71 2 1 80 90 1025 0
## 223 1 116 2 76 1 1 80 80 NA 0
## 224 1 188 1 77 1 1 80 60 NA 3
## 225 13 191 1 39 1 0 90 90 2350 -5
## 226 32 105 1 75 2 2 60 70 1025 5
## 227 6 174 1 66 1 1 90 100 1075 1
## 228 22 177 1 58 2 1 80 90 1060 0Fitting the model
survival_function = survfit(Surv(lung$time, lung$status == 2) ~ 1)
survival_function## Call: survfit(formula = Surv(lung$time, lung$status == 2) ~ 1)
##
## n events median 0.95LCL 0.95UCL
## [1,] 228 165 310 285 363Plotting the function
plot(survival_function, main = "Kaplan-Meier", xlab = "Number of days",
ylab = "Prob of survival")
Cox proportional hazard model
Uses the hazard function, which considers independent variables in regression.
Also does not assume an underlying probability, but assumes that the hazards are constant over time.
More useful to measure instantaneous risk of deaths, often delivers better results because more volatile with data and features.
Tends to drop sharper as time increases.
cox_mod <- coxph(Surv(lung$time, lung$status == 2) ~., data = lung)
summary(cox_mod)## Call:
## coxph(formula = Surv(lung$time, lung$status == 2) ~ ., data = lung)
##
## n= 167, number of events= 120
## (61 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## inst -3.037e-02 9.701e-01 1.312e-02 -2.315 0.020619 *
## age 1.281e-02 1.013e+00 1.194e-02 1.073 0.283403
## sex -5.666e-01 5.674e-01 2.014e-01 -2.814 0.004890 **
## ph.ecog 9.074e-01 2.478e+00 2.386e-01 3.803 0.000143 ***
## ph.karno 2.658e-02 1.027e+00 1.163e-02 2.286 0.022231 *
## pat.karno -1.091e-02 9.891e-01 8.141e-03 -1.340 0.180160
## meal.cal 2.602e-06 1.000e+00 2.677e-04 0.010 0.992244
## wt.loss -1.671e-02 9.834e-01 7.911e-03 -2.112 0.034647 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## inst 0.9701 1.0308 0.9455 0.9954
## age 1.0129 0.9873 0.9895 1.0369
## sex 0.5674 1.7623 0.3824 0.8420
## ph.ecog 2.4778 0.4036 1.5523 3.9552
## ph.karno 1.0269 0.9738 1.0038 1.0506
## pat.karno 0.9891 1.0110 0.9735 1.0051
## meal.cal 1.0000 1.0000 0.9995 1.0005
## wt.loss 0.9834 1.0169 0.9683 0.9988
##
## Concordance= 0.648 (se = 0.03 )
## Likelihood ratio test= 33.7 on 8 df, p=5e-05
## Wald test = 31.72 on 8 df, p=1e-04
## Score (logrank) test = 32.51 on 8 df, p=8e-05Fitting the model and plotting
cox <- survfit(cox_mod)
plot(cox, main = "Cox proportional hazard model", xlab = "Number of days",
ylab = "Prob of survival")
More advanced example of survival models
Source: https://rviews.rstudio.com/2017/09/25/survival-analysis-with-r/
Loading the data and packages
library(survival)
library(ranger)
library(ggplot2)
library(dplyr)
library(ggfortify) ## for better survival plots
data(veteran)## Warning in data(veteran): data set 'veteran' not foundhead(veteran)## trt celltype time status karno diagtime age prior
## 1 1 squamous 72 1 60 7 69 0
## 2 1 squamous 411 1 70 5 64 10
## 3 1 squamous 228 1 60 3 38 0
## 4 1 squamous 126 1 60 9 63 10
## 5 1 squamous 118 1 70 11 65 10
## 6 1 squamous 10 1 20 5 49 0Kaplan Meier Analysis
km <- with(veteran, Surv(time, status))
head(km, 80)## [1] 72 411 228 126 118 10 82 110 314 100+ 42 8 144 25+ 11
## [16] 30 384 4 54 13 123+ 97+ 153 59 117 16 151 22 56 21
## [31] 18 139 20 31 52 287 18 51 122 27 54 7 63 392 10
## [46] 8 92 35 117 132 12 162 3 95 177 162 216 553 278 12
## [61] 260 200 156 182+ 143 105 103 250 100 999 112 87+ 231+ 242 991
## [76] 111 1 587 389 33Fitting model with estimates at 1, 30, 60, and 90 days
km_fit <- survfit(Surv(time, status) ~ 1, data = veteran)
summary(km_fit, times = c(1, 30, 60, 90*(1:10)))## Call: survfit(formula = Surv(time, status) ~ 1, data = veteran)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 137 2 0.985 0.0102 0.96552 1.0000
## 30 97 39 0.700 0.0392 0.62774 0.7816
## 60 73 22 0.538 0.0427 0.46070 0.6288
## 90 62 10 0.464 0.0428 0.38731 0.5560
## 180 27 30 0.222 0.0369 0.16066 0.3079
## 270 16 9 0.144 0.0319 0.09338 0.2223
## 360 10 6 0.090 0.0265 0.05061 0.1602
## 450 5 5 0.045 0.0194 0.01931 0.1049
## 540 4 1 0.036 0.0175 0.01389 0.0934
## 630 2 2 0.018 0.0126 0.00459 0.0707
## 720 2 0 0.018 0.0126 0.00459 0.0707
## 810 2 0 0.018 0.0126 0.00459 0.0707
## 900 2 0 0.018 0.0126 0.00459 0.0707Plotting
autoplot(km_fit)
Survival curves by treatment
km_trt_fit <- survfit(Surv(time, status) ~ trt, data = veteran)
autoplot(km_trt_fit)
Plotting patients by age
vet <- mutate(veteran, AG = ifelse((age < 60), "LT60", "OV60"),
AG = factor(AG),
trt = factor(trt, labels = c("standard", "test")),
prior = factor(prior, labels = c("NO", "Yes")))
km_AG_fit <- survfit(Surv(time, status) ~ AG, data = vet)
autoplot(km_AG_fit)
Cox analysis
Fitting model
cox <- coxph(Surv(time, status) ~ trt + celltype + karno + diagtime + age
+ prior, data = vet)
summary(cox)## Call:
## coxph(formula = Surv(time, status) ~ trt + celltype + karno +
## diagtime + age + prior, data = vet)
##
## n= 137, number of events= 128
##
## coef exp(coef) se(coef) z Pr(>|z|)
## trttest 2.946e-01 1.343e+00 2.075e-01 1.419 0.15577
## celltypesmallcell 8.616e-01 2.367e+00 2.753e-01 3.130 0.00175 **
## celltypeadeno 1.196e+00 3.307e+00 3.009e-01 3.975 7.05e-05 ***
## celltypelarge 4.013e-01 1.494e+00 2.827e-01 1.420 0.15574
## karno -3.282e-02 9.677e-01 5.508e-03 -5.958 2.55e-09 ***
## diagtime 8.132e-05 1.000e+00 9.136e-03 0.009 0.99290
## age -8.706e-03 9.913e-01 9.300e-03 -0.936 0.34920
## priorYes 7.159e-02 1.074e+00 2.323e-01 0.308 0.75794
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## trttest 1.3426 0.7448 0.8939 2.0166
## celltypesmallcell 2.3669 0.4225 1.3799 4.0597
## celltypeadeno 3.3071 0.3024 1.8336 5.9647
## celltypelarge 1.4938 0.6695 0.8583 2.5996
## karno 0.9677 1.0334 0.9573 0.9782
## diagtime 1.0001 0.9999 0.9823 1.0182
## age 0.9913 1.0087 0.9734 1.0096
## priorYes 1.0742 0.9309 0.6813 1.6937
##
## Concordance= 0.736 (se = 0.021 )
## Likelihood ratio test= 62.1 on 8 df, p=2e-10
## Wald test = 62.37 on 8 df, p=2e-10
## Score (logrank) test = 66.74 on 8 df, p=2e-11Plotting
cox_fit <- survfit(cox)
autoplot(cox_fit)
Notes on cox regression models
Need to be cautious when interpreting results because it assumes covariates are constant, and in many cases there are not
Hazard: the rate at which events happen. CPH assumes hazards are constant over time
Hazard ratio: ratio of hazard between covariates (such as male/female)
Look more into Concordance stats if wanting to use ROC curves to assess model performance
Aalen’s additive regression model shows how covariates change over time
aa_fit <- aareg(Surv(time, status) ~ trt + celltype +
karno + diagtime + age + prior, data = vet)
aa_fit## Call:
## aareg(formula = Surv(time, status) ~ trt + celltype + karno +
## diagtime + age + prior, data = vet)
##
## n= 137
## 75 out of 97 unique event times used
##
## slope coef se(coef) z p
## Intercept 0.083400 3.81e-02 1.09e-02 3.490 4.79e-04
## trttest 0.006730 2.49e-03 2.58e-03 0.967 3.34e-01
## celltypesmallcell 0.015000 7.30e-03 3.38e-03 2.160 3.09e-02
## celltypeadeno 0.018400 1.03e-02 4.20e-03 2.450 1.42e-02
## celltypelarge -0.001090 -6.21e-04 2.71e-03 -0.229 8.19e-01
## karno -0.001180 -4.37e-04 8.77e-05 -4.980 6.28e-07
## diagtime -0.000243 -4.92e-05 1.64e-04 -0.300 7.65e-01
## age -0.000246 -6.27e-05 1.28e-04 -0.491 6.23e-01
## priorYes 0.003300 1.54e-03 2.86e-03 0.539 5.90e-01
##
## Chisq=41.62 on 8 df, p=1.6e-06; test weights=aalenautoplot(aa_fit)
Tuning cox proportional hazards regression
Source: https://www.youtube.com/watch?v=TrS2M5imOt8&t=9s
exp(coef)shows the hazard ratio- for quantitative variables, the hazard ratio compares a 1 unit difference
- Concordance: goodness of fit for coxph models. Percentage of cases that model gets correct (50% is random guessing)
- Try to improve concordance score by adding variables, testing nested models?
additional research on JMbayes2 model
Source: https://drizopoulos.github.io/JMbayes2/articles/JMbayes2.html
Primary function is jm()
Basic use example
Goal: to assess the strength of the association between the risk of death and levels of serum bilirubin
library(JMbayes2)
pbc2.id$status2 <- as.numeric(pbc2.id$status != 'alive')
CoxFit <- coxph(Surv(years, status2) ~ sex, data = pbc2.id)Describe patient profiles for biomarker using a linear mixed model
fml <- lme(log(serBilir) ~ year * sex, data = pbc2, random = ~ year | id) Linking survival and longitudinal models
jointFit1 <- jm(CoxFit, fml, time_var = "year")
summary(jointFit1)##
## Call:
## jm(Surv_object = CoxFit, Mixed_objects = fml, time_var = "year")
##
## Data Descriptives:
## Number of Groups: 312 Number of events: 169 (54.2%)
## Number of Observations:
## log(serBilir): 1945
##
## DIC WAIC LPML
## marginal 4363.693 5356.479 -3358.346
## conditional 3538.695 3357.098 -1901.875
##
## Random-effects covariance matrix:
##
## StdDev Corr
## (Intr) 1.0032 (Intr)
## year 0.1839 0.4006
##
## Survival Outcome:
## Mean StDev 2.5% 97.5% P Rhat
## sexfemale -0.1529 0.2694 -0.6499 0.3888 0.5636 1.0007
## value(log(serBilir)) 1.2507 0.0866 1.0810 1.4238 0.0000 1.0117
##
## Longitudinal Outcome: log(serBilir) (family = gaussian, link = identity)
## Mean StDev 2.5% 97.5% P Rhat
## (Intercept) 0.7238 0.1720 0.3852 1.0617 0.0000 0.9998
## year 0.2673 0.0382 0.1930 0.3445 0.0000 1.0037
## sexfemale -0.2639 0.1824 -0.6187 0.0887 0.1527 0.9998
## year:sexfemale -0.0887 0.0407 -0.1687 -0.0094 0.0276 1.0046
## sigma 0.3465 0.0065 0.3342 0.3597 0.0000 1.0105
##
## MCMC summary:
## chains: 3
## iterations per chain: 3500
## burn-in per chain: 500
## thinning: 1
## time: 7 sec- Assume the instantaneous risk of an event at a specific time \(t\) is associated with the value of the linear predictor of the longitudinal outcome at the same point \(t\)
Traceplots:
MCMC stands for Markov chain Monte Carlo
Traceplots are used to determine how many iteration it will take to summarize the posterior distribution (convergence) (think Bayes and conditional probabilities)
- If model converges, then traceplot moves around the mode of distribution
ggtraceplot(jointFit1, "alphas")
Density plot
ggdensityplot(jointFit1, "alphas")